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Abstract 

Epidemic dynamics in a stochastic network of interacting epidemic centers is considered. 
The epidemic and migration processes are modelled by Markov’s chains. Explicit formulas 
for probability distribution of the migration process are derived. Dependence of outbreak 
parameters on initial parameters, population, coupling parameters is examined analytically 
and numerically. The mean field approximation for a general migration process is derived. An 
approximate method allowing computation of statistical moments for networks with highly 
populated centres is proposed and tested numerically. 
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Introduction 

Epidemic outbreak in a populated center or in a network of populated centra develops stochasti¬ 
cally due to random interaction between individuals inside a center and due to a random migration 
between centers of the network. Conventionally, an epidemic outbreak in a highly populated cen¬ 
ter is described by deterministic processes in according with Law of Large Numbers (LLN) [[Q. 
The mean field approximation (hydrodynamic limits) of the appropriate statistical models estab¬ 
lishes the basic relation of stochastic description to the dynamical equations, say the SIR model 
(susceptible/infected/removed) and its more sophisticated modifications (SEIR, SIS, MSIR, etc.) 

However, there are two important cases when stochastic effects are essential. Firstly, it is ob¬ 
viously important when the populations in centers are not large. The second less obvious scenario 
can occur at the initial stage of outbreak when the number of infectives is small, then the discrete¬ 
ness of population can essentially affect the dynamics of the outbreak making it stochastic. For an 
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isolated center, this case is thoroughly studied in |0. Analysis of a network of interacting epidemic 
centers requires an account of migration fluxes between them. 

So, if the initial number of infectives triggering the outbreak in a particular populated center is 
small (that is typical for many outbreaks) than the LLN fails at least at initial stage until the number 
of infectives is large enough. For this reason the observed number of infectives can be significantly 
different from the prediction of a deterministic model, i.e. the standard deviation of the number of 
infectives and the peak outbreak time can be wide even in highly populated center or a network of 
such centers. 

In principle, the probability density function (PDF), its standard deviation, and other important 
characteristics for the outbreak forecast could be determined by a direct numerical simulation. 
However, this simulation may be very costly and require too much CPU time even for modem 
computers. Our goal is to develop a technique for a rough estimation of the outbreak statistical 
characteristics skipping such huge computation by applying some perturbation methods. 

Our toolkit is the so-called small initial contagion (SIC) approximation for the case of a large 
populated center when initial number of infectives is small, cf. |Z| in the case of a single populated 
center. For a network of highly populated SIR centers in the framework of deterministic model 
such technique is described in 0141. In this paper we develop a stochastic version of SIC approach 
based on the assumption that in the real epidemic centra the number of infectives triggering an 
outbreak is still small. In a sense, the SIC approach looks like a key to solving cumbersome 
epidemic networks. 

In 0, a randomized analogue of the standard SIR model is considered. In the SIC approxima¬ 
tion, one distinguishes two linked stages of epidemic evolution. At the stage 1 of initial contami¬ 
nation the number of infectives is small and a discrete formulation is vital. At this stage the system 
is randomized, and governed by stochastic equations. At the stage 2 of developed outbreak when 
the numbers of individuals in all the components are large the LLN works and the standard deter¬ 
ministic SIR model can describe the outbreak process accurately enough. Therefore it is natural to 
consider a deterministic system with randomized initial conditions linked to the stochastic stage 1. 
The statistical characteristics of the complete model are obtained by applying deterministic equa¬ 
tions with random initial conditions using the matching asymptotic expansions technique (cf. If5l). 
In contrast to the traditional technique, the asymptotic approximation of a randomized evolution (at 
a brief initial period) is matched with a deterministic evolution with randomized initial conditions 
(for all other times). Nevertheless, as in the traditional approach, we match the approximations at 
some intermediate time t in the interval where the both approximations are valid and which drops 
out from the final results (cf. 0). 

The Markov chain (MC) describing the randomized SIR model has been studied previously, 
e.g., in [j6l where a partial differential equation (PDE) for the moment generating function was de¬ 
rived. We develop a similar approach for a pair of linked centers and obtain approximate formulae 
for their main statistical characteristics. The results of large-scale numerical simulation are in a 
good agreement with the appropriate models. 

The paper is organized as follows. In Section Q] we introduce a MC model describing random¬ 
ized epidemic outbreak in two populated centers coupled by a random migration of all types of 
species. Here we consider convergence of this model to a deterministic mean-field model pro- 


2 


posed in 0. In Section [2] we describe a model of random migration between two interacting 
SIR centers taking place before the outbreak started to determine all the initial conditions for the 
MC model. Migration between centers is also described as a Markov chain. Here we derive the 
Master/Kolmogorov’s equations for the probability generating functions (PGF) and solve them an¬ 
alytically. This analysis confirms the diffusion-like model of migration heuristically proposed in 
If7l . In Section [3] the numerical algorithm for direct solving the MC model is described, the de¬ 
pendence of outbreak parameters on the population size, initial number of infectives and migration 
parameters are presented and discussed. In Section |4] the generalization of the MC model on an 
arbitrary network of the epidemic centers is studied, and difficulties of direct numerical modelling 
are estimated. Here the simplified model which looks relevant for a network of a highly popu¬ 
lated centers is proposed. In Discussions we make a comparison with some previously considered 
models and outline the perspective of the future development. 


1 Randomized model of two SIR centra interaction 

Consider two populated nodes, 1 and 2, with populations Ni and N 2 , respectively. Let S n (t), 
I n (t), R n (t ) be the numbers of host susceptibles, infectives and removed, respectively, in node n 
at time t. Let S mn (t), I mn (t), R mn (t) be numbers of guest susceptibles, infectives and removed, 
respectively, in node n migrated from node m at time t. Note that in the standard SIR model, 
removed individuals do not interact with other species, do not affect dynamics of susceptibles and 
infectives, and can be omitted from consideration [jHHHU. 

Assume that the populations in every node is completely mixed, the contamination rate /3 n of 
a susceptible individual in node n at time interval \t, t + dt] is proportional to the number of all 
infectives in node n: host infectives I n at time t plus guest infectives I mn immigrated from node 
m. Next, every infective in node n can be removed with probability rate a n (cf. 0). 

The model is described by the migration rate 7 nm from node n to node m and return rate 5 mn 
for a guest individual to return to his host node, they may be different for different species, i.e. we 
specify the migration process for susceptibles by parameters 7 nm^mn an d for infectives— ll vm ^l nn 
(cf. 0). Obviously, return rate of the host node should be higher than the migration rate to a 
neighbour node, i.e., Y mn < 5' inV 7 ;L < Cn- 

Taking into account all above, we model a network of two SIR centers interacting due to migra¬ 
tion of individuals between them a Markov chain (MC) which full description is given in Table 1. 

If Jo infectives appear in center 1 at time t = 0 than the initial conditions for this process are 

h = Iq, 5', = TV, - 7 0 - S° 12 , I 12 = 0, S 12 = S? 2 , 

h = 0, S 2 = N 2 -S° 21 , I 2l = 0, S 21 = S° 2V { ) 

Here the initial numbers of guest susceptibles S® 2 and .S'^ are random and determined by mi¬ 

gration processes between centers taking place before appearing of a single infective. In Section [2] 
we determine this distribution (which turns out to be binomial) by considering pure migration pro¬ 
cesses taking place before the outbreak: see Eq. (fl9l) below. Mean values for S® 2 and are given 
by®. 
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Numerical simulation based on this model is presented and discussed in Section [3] Analyt¬ 
ical approach via Master/Kolmogorov equations (see j2j) is too cumbersome, their analysis and 
solution in the general case is rather complicated, so we will develop reasonable approximations. 


Table 1: Markov’s chain for two coupled SIR nodes 


i 

Process #?’ 

Rate, 

Restriction 

Description 

1 

Co 

■1 1 
Co 

+ 1 

I— 1 I— 1 

/3i(7i+7 2 i)Si 

7, < A, 

Contamination in 1 (host) 

2 

S21 £ S 2 1 — 1 
h\ ->■ / 2 i + 1 

A(7i+7 2 i)S 2 i 

72i < (V 2 

Contamination in 1 (guest) 

3 

/r^/r-1 

«i7i 


Recovering in 1 (host) 

4 

I'l\ —£ 7 2 l — 1 

a 1 7 2 1 


Recovering in 1 (guest) 

5 

5*! —>■ 5i -1 

7712 —£ Si2 + 1 

7i 5 2 Si 

S 12 < N x 

Migration from 1 to 2 

6 

h ~+h -1 
h-2 —■> I\2 + 1 

73 nh 

7ia 7 S7 

Migration from 1 to 2 

7 

Si ->S X +1 
Sl2 —> S 12 — 1 

4s 12 

Si < TVi 

Return from 2 to 1 

8 

Jl ->-7, +1 
h-2 —■> Il2 — 1 

7 2 17l2 

7, < A, 

Return from 2 to 1 

9 

l i 

+ 1 

h-* 

P 2 (7 2 T7i 2 ) S 2 

7 2 7 iV 2 

Contamination in 2 (host) 

10 

Si 2 —> s 12 — 1 

112 ~> 112 +1 

/3 2 (72+7i2)Si2 

7i2 7 iVi 

Contamination in 2 (guest) 

11 

I 2 £ 7 2 1 

ot 2 I 2 


Recovering in 2 (host) 

12 

I\2 ~> 112 ~ 1 

ot 2 h 2 


Recovering in 2 (guest) 

13 

S 2 —£ S 2 — 1 
S 2 I ““£ S 21 + 1 

T 21 S 2 

S 21 7 iV 2 

Migration from 2 to 1 

14 

7 2 “£ 7 2 — 1 
7 2 i —£ 7 2 i + 1 

721 J 2 

7 2 i 7 A^ 2 

Migration from 2 to 1 

15 

S 2 —£ S 2 +1 
S 21 —£ S 21 — 1 

5f 2 S 2 i 

s 2 <n 2 

Return from 1 to 2 

16 

h —£ 12 +1 

I 21 —> 7 2 i; — 1 

^V2^2\ 

I 2 7 £V 2 

Return from 1 to 2 


Proposition 1. The scaled Markov chain (MC) {/*(£), S*(t), 7* m (f), S* m (t)} (n — 1, 2, m — 2,1) 

m = A-'/„(*). S-Jt) = A-'S„((), 

/„(() = A-■/_(!), S’ m (t) = A-'U*l 

in a populations of sizes AN n ,AN m related with process (7 n (£), S n (f), 7 rem (£), S nm (£)} defined 


4 
























in Table 1 and subjected to initial conditions © by scaling the transition rates /?„—>• A 1 /3 n , and 
scaling of initial conditions as 

/ 1 (0) = A/ 0 , S 1 (0)=A{N 1 -I o -S° 2 ), 

/ 2 (0) = 0, S 2 (0) = A(N 2 -S° 21 ), 

/i 2 (0) = 0, S 12 = AS° 2 , 

/ 2 i(0) = 0, S 21 = AS° 21 


where the PDFs of independent random variables S® 2 and have binomial distributions (fl9l) . 
converges in distribution as A —» oo to the deterministic functions | I n (t), S n (t). I nm (t). S n 
satisfying ODEs introduced in I© 

—finSn{In + Imn ) — Inm^n + ^mn^nm (4) 

A?) S n ( I n T Imn) C^nln 7/im-7 4" ^rnnInm (5) 

PnSmn(In 4" Imn) 4“ Tmn^i AimS'mn (6) 


d c 
d t° n 


d t 


In = 


_d n _ 

a 4 - ^m,n 


d t 


An S nl n (In A Imn) O imn T 4m//. In 


Am A ran 


(7) 


subjected to the initial conditions 


A(o) = Jo, Si (0) — N { — I 0 — Sy 2l 

7 2 (0) = 0, S 2 (0) = N 2 -S° 21 , 

4 (o) = o, Si2 = sj 2 , 

I 2 i (0) = 0, S 2 \ = S 21 


(8) 


where 

c0 _ ^ 12^1 c»0 _ 7^2 

12_ 7f 2 + 4’ 21 “7 2 S ! + 5f 2 


(9) 


In fact, equations ©-© can be derived phenomenologically: if the number of species is large 
enough, its change by one or by few can be considered as infinitesimally small. For example, 
the number of infectives I\ can increase due to process #1 and #8 with rates Ai (I\+I>\)S\ and 
5 ,21 /1 2 , respectively, or decrease due to process #3 and #6 with rates a \ I \ and 7(3/1, respectively. 
Therefore the rate of dii in time interval d t can be estimated as d/i = \ji\ (7, +/ 2 | ),Sj + S 21 Ii 2 ] (it — 
[ai/i + 7( 2 /i] (it. that gives Eq. ([5]) for n — 1, m — 2. The same holds for all other species. This 
argument can be made rigorous with the help of Law of Large Numbers (LLN). Then it indicates 
that the mean values of the variables converge to the mean-field (hydrodynamic) limit. 

The more delicate question is to establish the convergence in probability. Sketch of the proof 
is given in Appendix. 

2 Random migration of non-contaminating species 

In order to elaborate distributions for S® m taking place at the very beginning of the outbreak we 
study the pure migration process setting I\ = 0, 1 2 = 0. In this case the MC described in Table 1 
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can be split into two independent processes: Si S \ 2 — N\ — S \, and S 2 +> S 2 \ = N 2 — S 2 . For 
each of them we have the following MC in terms of a single random variable S n , n = 1 , 2 : 


Process 

Rate 

Sn ^ S n 1 
Sn ^ S n + 1 

'V s ' <7 

Inm^n 

~ S„ ) 


( 10 ) 


Let Pk(t) = P(S nm (t)=k) = P (S n (t)=N n — k ) be the probability distribution in node m at 
instant t. Then Masters/Kolmogorov’s equations take the form 

j t P k = 7 (jN - k + 1]^) P k .i - 7 ( 1 V - k)P k + 5 (jk + 1]^) P k+ 1 - 5kP k (11) 

where 0 < k < N; for the sake of simplicity we temporary set 7 = 7 nm> 5 = N = N n • Here 
the following notation is introduced to write equations for k = 0 and k = N in the same form as 
others 


r ktf = 


k , 0 < k < N 

0 , otherwise. 


Equations (fill) implies the following PDE 

G t = {l-z) [(-7 z-S) Gz + 'yNG ] 
for the probability generating function (PGF) 


( 12 ) 


(13) 


N 


G(z,t) = J2 zkp k(t)- 


k =0 


The initial condition P o (0) = 1, P fe >i(0) = 0 implies 

G(z, 0) = 1. 

The solution to problem (fT3l)-(fl5l) can be found explicitly 

(7 z + 5 ) — 7 (z — l)e _ ( 7+<5 l t ’ 1 N 


G(z,t) = 


7 + 5 


Now one can easily calculate all the moments of distribution {P k (t)}, say 

ES{t) = m (t) = G z { 1 , t) = Ne [l — e~ t/T ] 
var (S(t)) = fi 2 (t) = G zz (l,t) + jii - n\ = /ii(t) [ee~ t/r + (1 - e)] 


(14) 


(15) 


(16) 


(17) 

(18) 


where e = 7/(7 + 5 ), r = 1/(7 + 5 ). 
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If the migration process lasts long enough before the outbreak starts then the PGF takes its 
limiting form for £ —>■ oo 

G(z, oo) = (ez + (1 — e)) N 

which is the MGF for a binomial distribution: 


N r . 


P(S“ m = k) = ( ) (4J (1 - 7,J 


S \Nn-k 


From here on we return to the indexed notation 


n,S,i 

gS,I _ inm 


S,I , c5,/' 

jnm i Omn 


This distribution has the following first two moments 


S nm = E *5 nm = /ri(oo) = e nm N n 
w(i) = ^ 2 ( 00 ) = N n e s nm (l - £ S nm ). 


(19) 


( 20 ) 


( 21 ) 

( 22 ) 


The relative standard deviation (i.e. for the process X = S nm /S® m ) decays as N n l/2 : 


CTcO /CO 
> ~ J nm / *-'n- 


^ 2 ( 00 ) 

^ 1 ( 00 ) 


(1 -e s ) 

V 0 nm ) 


e s N 


(23) 


Hence, when iV n —» 00 , the migration process tends in probability to the deterministic limit de¬ 
scribed in 0 . 

Thus, in the MC model defined by Table Q] the initial conditions S { \ 2 and S 2l can be selected 
randomly with the binomial distribution (fl9l) or approximated by a Gaussian function if N n is large 
enough. 

Parameter e s nm := 7 f m /( 7 + ^nm) = $nm/N n indicates the mean share of individuals from 
node n migrated to node rn. Obviously this share in average should be small for highly populated 
centers: say, half population of a city hardly can be on a visit into another city for an essential time. 
Parameter £ s nrn can be treated as a coupling parameter characterizing how intensive are migration 
fluxes between populated centers. Analogous fluxes of infectives hardly are more intensive, there¬ 
fore = 7 2 m /( 7 2 m + 5 T nm ) < e s nm . So, for highly populated centers, the following inequality 
should be kept 

4m < 1 ^ 7mn < 4m- ( 24 ) 

The second important characteristic of migration process is the characteristic migration time 
T nm = 1 /( 7 nm + 4 m)- F rom (Et one can see that it indicates how soon the dynamic equilibrium 
established after the migration process started or the population changes suddenly. Both pairs of 
parameters: { 7 , 5} and (e, r} are uniquely related. 
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3 Direct numerical simulation of two interacting SIR centra 


3.1 Numerical scheme 

In the numerical simulation of the randomized SIR model the time interval was divided into small 
steps At such that the sum of all rates from Table 1 multiplied by At is essentially less than 1: 


max {r's(t)} At <C 1 ==>• At = min {Pt/v^(t)} 


(25) 


where P t is the admitted threshold, say, P t = 0.1. 

Probability that at least one event occurs in unit of time is bounded by the sum of rates of all 
the processes z/ s (t) = Y^=i u i '■ 

Vs(t) — Pi (h + hi) (Si + S 21 ) + P 2 (h + ^ 12 ) (S 2 + S 12 ) 

+ ot 2 (h + I 12 ) + oti(h + hi) 

+ 7^1 + 7l2-^l + $2lSl2 + ^ 21-^12 + I 21 S 2 + l2lh + $2lSl2 + £>l 2 hl 
In this relation we majorize h, Si < N u I 2 , S 2 < N 2 , h 2 , S 12 < N u I 2 1 , .5 2 1 < N 2 , then 

max(r's) — (Pi + /?2 ) (-A7i + N 2 )~ + (on + q? 2 )(Ai + N 2 ) 

+ (712 + 7l2 + ^21 + $2l)Nl + (721 + 721 + ^12 + ^2)^2 
+ 7^^! + 7lVl + $2lSl 2 + & 21 h 2 

Actually this value overestimates the really occurred total rate significantly as it is very improbable 
that numbers of guest susceptibles and infectives in a highly populated center essentially exceeds 
values £nm^n and £ T nm N n , respectively, where sS T n is defined in (fl9l) . in virtue of (1241) . For this 
reason we can account that S nm < e^ m N n and I nm < £^ m N n (and also use the rigorous inequalities 


I n S n < jN%). Then we obtain the more realistic estimation: 

max (u-s) ~ \PiNl + |/3 2 1V 2 2 


+ Pi ( £ 2 i + £ 2 i) N2N-1 + p 2 (e ( 2 + ef 2 ) ^2 Ni 
+ OL\ (Ni + £ 21 ^ 2 ) + a 2 (-A 2 + £ 12^l) 

+ (t?2 + 712 ) + (efi^ + £ 21 ^ 12 ) Ni 

+ (721 + 721) N 2 + (^2^21 + £12^21) N 2 

+ Pl £ 21 £ 21 N 2 N 1 + P 2 £ 12 ^ 12 -^2 A 1 

The following numerical scheme is used: 

1. Assign the initial values to 8 variables 

h = I 0 , Si = Ni-I 0 -S° 12 , I 12 = 0, S 12 = S° 12 , 
h = 0, S 2 — N 2 —S 21 , I 2 i = 0, S 21 = 

where Si 2 , S 21 are random numbers distributed in accordance with (fl9l) . 
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2. Calculate the current rates {z/j, i — 1,..., 16} indicated in Table 1. 

3. Calculate the current probability of at least one event occurrence in accordance with eq. (1251) : 

16 

p = A tvY.it) = At Vj 

i= 1 

4. Generate uniformly distributed random number x € [0,1]. 

5. If x > p then no events occur. In this case: 

(a) advance at one time step: t 4— t + At without changing variables h, ■ ■ ■, S 2 1 , also 
z'i, • • •, Cl6 and p. 

(b) if t > t rnax terminate the process, otherwise go to step 4. 

If x < p then at least one event occurs. In this case: 

(a) calculate the intervals Ay t = [r/,_,, ry], ry = ]T” =1 z/ ? 

(b) generate the second random number y uniformly distributed in [0, ^ 16 ] 

(c) find in which interval y falls. 

(d) perform the process described in Table 1 with the correspondent rate 

(e) advance at one time step: t 4—t + At 

(f) if t > f max terminate the process, otherwise go to step 2. 

3.2 Numerical results 

For numerical computation a basic model with two identical centers has parameters: A" 1j2 = 10 4 , 
qli j2 = 1, Roi ,2 = 4, £[’2 = 0.01, t [2 = 5 where Ro 1)2 = (^i, 2 /cii, 2 ) A 1>2 are reproduction 
numbers. The initial number of infectives in the first node J 0 = ioNi where i 0 = 0.01 is taken for 
the basic model. 

A few realizations of numerical computation are depicted in Figure |T] Here the time variations 
of the total number of infectives 1^ — I n + I mn in every node are shown and compared with the 
curves based on integration of the deterministic initial value problem dU)-© 

In the first set of numerical experiments the total population size varies from A r i = /V 2 = A r = 
400 up to 10 6 . Number of realization L was taken 10 4 mainly but the number of realization is 
taken greater for small population N = 400 and 2000 and smaller for extremely high populations: 
250k and 1000k. The current mean number of total infectives /j} 2 (t) and standard deviation are 
computed. The results of the first set are shown in Figure [2] Observe that the mean value of the 
random process (thin lines) converges to the solution of correspondent deterministic problem (bold 
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Time 

Figure 1: Examples of realizations of two randomized SIR models. The total number of infectives 
is plotted in node 1 (if — I\ + I 2 i ) by thin grey lines and in node 2 (If — I 2 + I 12 ) by thin black 
lines. Bold dashed lines indicates the hydrodynamic limit. (Ni = N 2 = 2 k, R 01 = Ro 2 = 4, 
e — 0.01, r = 5, Iq/Ni = 0.01) 




Figure 2: Evolution of the mean values for lf n /N m m = 1,2 (left) and their standard deviations 
(right). Grey curves for node 1, black lines for node 2. Dashed lines indicate the hydrodynamic 
limit described by eqs. ©-([ 7 ]). The node population is indicated near the top of the correspondent 
curve. 

lines). But the convergence is much slower for node 2: for N\ = 10k the mean trajectory prac¬ 
tically coincides with the deterministic limit, on the other side the same effect in node 2 requires 
N 2 ~ 250k. 

The convergence rate is examined in Figure [3} One can see that for node 1 the convergence 
rate almost coincides with 0(N~ 1 ^ 2 ), as for node 2 the decay rate tends to 0(N~ 1 ^ 2 ) only for 
sufficiently large population: N = 10 6 . Thus for the secondary contaminated node the account of 
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I s 

randomness is essential even if its population is large provided that the migration parameters e 1 ’ 2 
are small (0.01 in this case). Say, if iV 2 = 400 the standard deviation exceeds the mean value up to 
the outbreak time. 



Figure 3: Maximal value of the standard deviation for processes I\/Ni and I 2 /N 2 vs population 
N — JVj — N 2 . Black curves for node 1, grey lines for node 2. The slope of dashed line corre¬ 
sponds to the decay law A r ~ 1/2 . 


In the second set of numerical experiments we study dependence of mean number of infectives 
and their standard deviation from the initial number of infectives J 0 in node 1. It varies from 1 to 
100 (the share i 0 = Iq/Ni varies from 10 -4 to 1CT 2 ). The results are plotted in Figured Because 




Figure 4: Left: Mean values for if/Ni (black) and its standard deviation stdlf/Ni (grey) for 
different J 0 . The initial number of infectives in node 1 is indicated near top of the corresponding 
curve. Dashed line indicates the deterministic limiting solution. Right: the same for I 2 /N 2 . 

the outbreak time depends on initial number of infectives, the mean-field curves become essentially 
different. Therefore it is appropriate to shift the time so that the peak outbreak for different initial 
condition is at the same instant, say, t = 0. Then all the curve are very close to each other and 
practically coincide with curve for the limiting solution introduces in (3l|4j. Observe that the 
smaller the number of initial infectives the greater is the standard deviation (std) and the larger 
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is the difference between the mean curve for randomized process and the mean-field curve. Also 
observe that for node 1, the discrepancy of mean number of infectives from the hydrodynamic 
limit as well as the standard deviation monotonically decay with the growth of J 0 . In node 2 the 
analogous discrepancy and the std slightly change when the number of initial infective varies from 
10 to 100. 



Figure 5: Maximal value of the standard deviation for I\ /N\ and I 2 /N 2 process vs population J 0 
Grey curve is for node 1, black curve is for node 2. The dashed line has the slope corresponding 
to the decay law iV -1 / 2 . 


In the third set of numerical experiments we study dependence of the mean number of infectives 
on the coupling coefficient e (the same for all species). It was varied in the range 10 -4 ,..., 10 _1 . 
The results are plotted in Figures [6] and [7] Observe that e practically does not affect the standard 
deviation of total number of infectives in node 1. Discrepancy of the mean curve from the mean- 
field curve becomes noticeable only for small e\ e < 0.05. 



0.15 

C 01 

S 0.05 


-1 


0 

t — t. 
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Figure 6: Mean values for If/Ni (left) and its standard deviations std/f/A r i (right) for different 
migration coefficient e. Dashed lines indicate the mean-field limit. The initial number of infectives 
in node 1 is indicated near top of the corresponding curve. 
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Figure 7: Mean values for lf/N 2 (left) and its standard deviations std lf/N 2 (right) for different 
migration coefficient s. Dashed lines indicate the mean-field limit. The initial number of infectives 
in node 1 is indicated near top of the corresponding curve. 


As for the second node, the standard deviation grows monotonically with the decrease of the 
coupling. That indicates the importance to account for randomness of the epidemic process in the 
case of weak coupling (i.e. in the case of relatively slow migration fluxes). 



Figure 8: Left: Maximal value of the standard deviation for processes I\/N\ and I 2 /N 2 vs pop¬ 
ulation N = Ni = N 2 . Grey curve is for node 1, black curve is for node 2. The dashed line 
has the slope corresponding to the decay law N~ l I 2 . Right: Outbreak value for I\ /N\ and I 2 /N 2 
processes vs population N = N\ = N 2 Grey curve is for node 1, black curve is for node 2. The 
dashed line are for the mean-field values. 
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4 Two-stage semi-randomized model 

The MC model for two coupled SIR centra can be readily generalized for an arbitrary network of 
M mutually interacting SIR centers: 


Process 

S n —> S n — 1, I n —> / n +1 

*S mn ^ ‘S'inn 1 ■ I-mn ^ Imn Tl 

In —■ y In~ 1 

Ini n ^ Irn ri 1 

*S'n ^ S'/ / 1, S nrn )■ S nm T1 

In ^ In 1 1 Inm t In ri i T1 

Sn —> S'n+l ■) S nm ^ S nm 1 

T T i-1 T r -1 


Pn(In~^~ I mn ) ‘S'n 

Ai(4i + Yl m Imn)S„ 

T 


O^nln 

Glnlnin 

Q 

lnm u n 

ry 1 / 

lnm J -n 


where m,n — 1, M, m ^ n. 

Generically, to study migration fluxes (7 > 0) between every pair of centers (although 
of different rates) we have to consider 2 M host and 2 M(M — 1) guest species. The correspondent 
MC model will contain 4 M 2 fluxes. The total rate can be evaluated via z/v = 0(M 2 N 2 ). Therefore 
for a network containing significant number of highly populated centra (say, main cities in the UK), 
the time interval At will have to be taken extremely small, and, hence, the GPU time for a single 
realization will be too long to get the necessary statistics. Thus the MC should be simplified for 
numerical modelling. 


4.1 The small initial contagion (SIC) approximation 

The SIC approximation is based on the assumptions which look relevant for a network of highly 
populated centers: 

• Population in every center is high: N n 1. 

• Migration fluxes between the centers are small: <7 1. 

• Initial number of infectives in the first contaminated center (say, n = 1), is small: 7 0 <C A r i. 

• Reproduction number exceeds unity and is not close to it in all the nodes: Ro n := f3 n N n /a n > 
1 + r where r = 0(1), r > 0. 

In these assumptions, the outbreak process in every center can be split into the following main 
stages: 

1. Contaminating stage: number of infectives is small I n <7 N n , S n ~ N n and the fluxes of 
infectives caused by migration are essential for the outbreak process (except the first node). 
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2. Developed outbreak: I n 1, when the contribution of migration fluxes is negligible, (also 
the mean-field description for every individual realization looks adequately). 

3. Recovering stage: the node is not affected by infective immigrants and slightly affects con¬ 
tamination of other nodes. 

From these assumptions follows that the outbreak dynamics in the first node can be considered 
independently, and can be described by the following MC 


Process 

Rate 

S\ —> 5*i—1, I\ —y /i+l 

faSih 

h h-l 

afai 


(27) 


with the initial condition fi(O) = J 0 , Si = Ni — I 0 studied in JUiU. At the contamination stage 
we have Si ~ N l and the MC can be further simplified 


Process 

Rate 


l\ —y /1+1 

{faNi) h 

(28) 

h -)■ h-l 

atih 



Epidemic dynamics in node 2 at the contamination stage (S 2 ~ N 2 ) can be described by 
analogous MC with additional flux v{t) of infectives migrated from node 1: 


Process 

Rate 

I 2 —y / 2 +1 
I 2 —y I 2 —1 

(faN 2 )I 2 + v(t) 
Oi 2 I 2 


(29) 


with / 2 (0) = 0 where 

fat) = {fa N 2 )I 12 + S( 2 I 2 1 . (30) 

At the contamination stage processes /12(f) and hi{t) are practically independent of process 
I 2 (t). Thus we have to consider MC (l29l) with a random flux v{t) which statistics will be specified 
later. 


4.2 Calculation of moments 

First we consider a single realization of the flux u(t) and treat it as a deterministic function. Later 
we will use averaging on u{t) to calculate moments of distribution for I fat). 

4.2.1 Calculation of PGF for a single realization G(z, t). 

Let P k {t ) = P(/ 2 (f) = k) be the probability of k infectives I 2 at instant t. Initial condition is 

P 0 (0) = 1, P k > o(0) = 0. (31) 


15 




















Kolmogorov’s equations for MC (l29l) are 

±P k = P' 2 (\k - 11 - (J' 2 kP k + a 2 {\k + l]" 2 )P k+ i - a 2 kP k + uP k ._! - vP k . (32) 


Here [P 2 = (3 2 N 2 . For the PGF G(z, t ) tending iV 2 —> oo (with f3 2 —> 0, (3' 2 — const ) we obtain the 
following PDE 

Gt = (z — 1) [((J 2 z — a 2 ) G z + u(t)G\ (33) 

with the initial condition G(z, 0) = 1. Its solution can be written in the integral form 


G(z, t ) = exp< — 


X 2 (z— 1) v ( t ') d t! 


P 2 (z~l) ~ {(3' 2 z-a 2 ) eW-t) 


(34) 


where A2 = ft 2 — a 2 is the initial growth rate of infectives in the deterministic SIR model in the 
limit Iq/N —* 0 (limiting solution introduced in []3l 3J]). 


4.2.2 Calculation of first moment E [I 2 (t)]. 

The first conditional moment Hi(t \ v) = E (I 2 (t) \ u) for fixed u(t) is 


Hi(t | u) — G z (l, t) — f u(t')e X2 ^ 4 3d t' (35) 

Jo 

Averaging over all realizations for u(t) (with a time varying PDF f v {t )): jiiit) = E fii(t \ v) 





fv(t')is{t') e^-Gdt' 


du 



(36) 


where u(t) = E v{t). Thus the average number of infectives in node 2 in the contamination stage 
relates with the flux u(t) via the convolution 


E [I 2 (t)] = Hi (t) = i/(t)* e Mt . 


(37) 


4.2.3 Calculation of second moment var [I 2 (t)\. 

To calculate it we apply the Law of Total Variation (e.g. Ifljl): 

var [ I 2 (t )] = E [var (I 2 (t) | v)] + var [E (I 2 (t) \ v)]. 

1. The first addend in (l38l) can be found through the PGF G(z, t): 

var (I 2 (t) | v) = G zz (l,t) + Hi(t \ v) ~ \ v) 


2 ^ /“ 
■^2 Jo 


v(t') 


e 2A 2 (t-f) _ e A 2 (t-f) 


d t' + Hi(t | v)- 


After the averaging through v we obtain 


2^ 


E [var (I 2 (t) | v)\ = Hi(t) + / z'(f') 

Jo 


g2A _ e A 


dt 1 . 


(38) 
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( 39 ) 


Thus a first addend in (1381) can written as a sum of two convolutions 


E [var(/ 2 (f) | v)\ = ^ -is(t ) * e 2Az * + (1 — * e X2t 

A 2 A 2 

where 9{t) = P' 2 I l2 {t) + 8{ 2 I 2 i(t). 

2 . To calculate the second addend in (l38l) . we temporary add /x 2 to it. Now it can be expressed 
via the covariance of flux z/(t): 


var [E(/ 2 (f) 11/)] + p\{t) = E yJ v (0 e A2(t_t,) d t' 

= f f* E [v (t') v (t")] e A2(t_t,) df , e A2 ^ _t "Mf // . 

Jo Jo 

The function in the integrand can be represented as a sum 

E [u ( t') v it")] = cov [v ( t' ), v ( t ''')] + p(t')D{t"). 
Integration of the second addend in (l4ll) gives just the temporary added term 

Jo Jo 

Thus the second addend in (l38l) can be written through the following integral 

var [E(/ 2 (f) \v)]= [ [ cov [v (t 1 ), v it")] 


(40) 


(41) 


(42) 


o jo 


in which we have to calculate the covariance of flux u(t). 

If flux v{t) is a random process controlled by a MC, calculation of its covariance is a com¬ 
plicated task, and consideration of this needs a separate work. Remind that the flux is a linear 
combination of two MC processes (l30l) : v(t) = /3 2 /i 2 (f) + 8[ 2 hi{t)- Here we approximate / 12 (f) 
and/ 2 i(f) by two mutually independent Poisson processes with variable rates ^/ 12 (£) and it), 
respectively, where / 12 (f) and I 2 \ (t) are calculated below. In this approximation, using the inde¬ 
pendence of increments of the inhomogeneous Poisson flow (e.g. iflOl ) we can write 


cov [ 1 12 (*'), /12 (£")] ~ /i 2 (min {/', /"}), 

cov [1 2 1 (£'), I 2 i it")] ~ hi (min {£', t"}), 

COV [I\2 it '), hi {t")] » 0 . 


We justify this approximation numerically below. Thus, for the covariance of the flux we have 

cov \v it '), V (/")] = (min (f',£"}), xu(t) = i/3 2 f Iuit) + (8[ 2 ) 2 J 21 (f). (43) 

Equation (l43l) holds true because the second central moment of the Poisson process coincides with 
the first moment. In this approximation, it is sufficient to compute the first moment of flux v(t) in 
order to compute the second moment for J 2 . 
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With the account for (l43l) . we split integral (l42l) into two parts: 


var [E(/ 2 (£) | v)\ 


Jj J\ + J 2 

Ji = [ [ w(f /, )e A2(t - i " ) df"e A2(t - f,) df / 


'0 

r*t 


J 2 = 



zu(t')e X2 ^~ t "' > dt"e X2 ^ t ~ t ' df'. 


o a*' 


Integrating ,J\ by parts and ,J 2 directly we obtain that they both give the same answer 

J1 — J2 — * e 2Az< — -^-w(t) * e Azt . 

A 2 a 2 


Finally combining the above results we have 


var(J 2 ) = 


4 (^ 2 ; f 2A 2 t 


4(/3 : ' a2 


a 2 

2)2 a 


I\2 * e 2 + 


Ji 2 * e Azt 


A 2 


a 2 


J 2 1 * e 


2A 2 i 


« , 2 (ay ,7 ■ 

"I i °12 


Ao 


Ao 


J 2 i * e Azi . 


(44) 


4.2.4 Computation of the average flux u(t). 

It is natural to split the total flux into two parts v(t) = j 3 ’ 2 Ii 2 + S[ 2 I 2 1 = 172(f) + v 2 i{t). Flux 
process 172(f) described by the MC 


Process 

Rate 

/l2 —> /l2+l 

7i 2 /i 

Jl2 —> 112 1 

(^21 + « 2 ) I 12 


(45) 


(with /i 2 (0) = 0) coincides with that described by (l29l ) if we set /3 2 <— 0, u(t) 7- 7(3/1, a 2 7- 
(5( x + a2) • Then we can immediately write down a solution for the PGF 


pt 

G(z, t) = exp(7(2 (z- 1 ) h it') e -( <5 2i+ Q2 )d'-d di '} 

Jo 


and the first moment 

I12 = 7(2 [ h (f) e-^i +«2)d-t') d t / . 

Jo 

Flux process z/ 2 i is more complicated and can be described by the following MC 


Process 

Rate 

S 2 1 —> S21 + 1 

721^2 

521 —> S21 — 1 


/21 —» I 2 I + 1, S21 ~^ ^21 — 1 

/ 4 i Ti 5 * 2 1 

I 2 1 ~> I 2 I “ 1 

5(2^21 


(46) 


( 47 ) 
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with the initial conditions 5Vi (0) = 4^2, I 2 i (0) = 0. If we split the third event into two inde¬ 
pendent events L n —> I 2 i + 1 and S 2 1 —>■ — 1 with the same rate, we can split MC (1471) into 

two MCs. The first MC describes migration of host susceptives from node 2 to node 1 and their 
possible removal due to contamination: 


Process Rate 

S21 —> S21 + 1 72^2 

5*21 —> S21 — 1 (<5f 2 + /3i/j) 5*21 

It is independent of the second MC with the rates 


Process Rate 

I 21 —> I21 +1 

121 —f I21 — 1 (^12 + 07) -^21 


(48) 


(49) 


which describes migration of susceptibles to a neighbor node, their contamination there and return 
to the host node as infected species. 

First, we study the first MC: (1481) . In accordance with (fl9l) it has the binomial initial distribu¬ 
tion: 

ft(o) = (^ 2 ) (4)‘ (1 - • (so) 

The probabilities P k (t ) = P(S' 2 i(f) = k) of k guest susceptibles 5 ' 2 i at instant f satisfy Kol¬ 
mogorov’s equations for MC (l29l) 



u(P k _ i; *P k )+a (\k+l]" 2 )P k+1 -kP h 


(51) 


where v = ^ 21 ^ 2 , « = (4 + Pih) • System (l5ll) implies the following PDE for MGF G(z, t) = 

E£L 0 z k P k (t) 

G t = {z - 1) [~a(t)G z + vG]. 


N' 


The initial condition is 

N 

GM) = y>- 

k =0 

The initial value problem (I52l)-(l53l) admits the explicit solution 

G(z, f) — [1 + e(z — 1 )(p{t)] N exp | u(z — 1 )0(f) 

where 0(f) = exp |— f* c»(f , )df , |. From here we have 


£ fe (i - ey~ k = (1 -E + £zy\ 


rt d t' \ 
) 00 ') 1 


ES 2 i(t)=G z (l,t) = N 2 


£ 21 + 721 / df'/0(f') 


0(f). 


(52) 


(53) 


(54) 


(55) 
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Now for process (l49l) in analogy with processes (l29l) and (l45l) we can immediately write 

hi = A [ e -(^ +Q2 ) (t - t ' ) E [h {t ') 5 2 i( 0] df'- (56) 

Jo 

Neglecting the mutual dependence of processes S 21 (t) and I \(t) we approximate 

E[S 21 (t)h(t)}^ES 21 (t)h(t). (57) 

Below we show numerically that it is satisfactory approximation for our applications. 

4.3 The second stage 

Remind that equations (l36l) . (l42l) . (1431) . (l44l) . (l46l) . (l56l) . (1571) are valid at the contamination stage 
only (S 2 ~ N 2 ). They allow us to calculate the first and second moments for the number of in- 
fectives without modelling the random process directly. To evaluate the moments at the developed 
outbreak we use the same approach as in [EJ: by approximating the outbreak via the mean field 
solution for a single SIR node with random initial conditions. 

For this aim we define the intermediate time t * such that the number of infective is large enough 
to use the mean field solution but the number of infective still slightly deviate from N 2 : J 2 (t*) 1 

and N 2 — S 2 -C N 2 . 

Then we generate L times a random number A" lognormally distributed (to guarantee the pos¬ 
itiveness) with mean J 2 (f*) and variance var(J 2 (t*)) and integrate the classical SIR equations: 
^S 2 = ~/3 2 S 2 I 2 , ±I 2 = (JJ 2 S 2 - a 2 )I 2 with initial condition / 2 (f*) = A", 

S 2 (U) = -W-! [—/3' 2 N 2 exp {fi 2 (X - N 2 )}} /& 

where Wf~[-] is the A th branch of the Lambert function IfTTlI . 

Let us emphasize that in the classical SIR model, the numbers of susceptives and infectives are 
related as / = iV — S + In [S'/ (N — J 0 )] //T where f3' = (JN (cf. (6l|). Resolving this relation with 
respect to S we can write 


S = -W k [-/3'(N - J 0 ) exp {/?'(/ - N)}} //T 

where k = —1 for the growing part and k = 0 for the decaying part of the outbreak. If the outbreak 
is triggered by a infinitesimal number of infectives we can set S = —W k [—/3'N exp {/3'(I — N )}] / / 3 '. 

Also note that it is natural to approximate the solution to a standard SIR model by the limiting 
solution ( h/N — > 0) introduced in J3|. The limiting solution is independent of the initial condition, 
therefore it is not necessity to integrate the ODEs L times but only once. 

Thus the proposed two-stage model of a coupled randomized epidemic centers allows us to 
calculate its first moments much faster than the direct simulation summarized in Table [T] 
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Table 2: Markov’s chain for two coupled SIR nodes in the SIC approximation 


# 

Process 

Rate 

1 

Si Si-1 

h ->• Ji+1 

PihSi 

2 

h h~l 

aJi 

3 

I\ 2 —> +L2 + 1 

Yuh 

4 

I\ 2 —> I\2 — 1 

5 21 Il 2 

5 

1S21 ~^ <521+1 


6 

S'21 —> S 2 i~ 1 

(5f 2 + Pih) s 21 

7 

I'21 —> /21 + 1 

jhIiS 2 i 

8 

hi ^21 1 

Qii hi 

9 

S 2 —> S 2 —1 
h ^2 + 1 

fi 2 {I 2 + Il2)S 2 + 5(2^21 

10 

I 2 —> I 2 — 1 

Oi 2 I 2 


4.4 Comparison with numerical simulation 

To show the accuracy of the proposed model we compare the solution obtained by different ap¬ 
proaches. Again the basic model of Section [3] is used: N l — N 2 — 10k, /3[ 2 — 4, a^ 2 = 1, 
efp = 0.01, t ^2 = 5, Iq/Ni = O.Olbut also model with the smaller population N\ = N 2 = 2k. 
We compare (i) the full randomized model described in Table [Tj which we regard as a benchmark; 
(ii) the SIC approximate randomized model where only flux of infectives from node 1 to node 2 is 
accounted, it is described in Table [2l (iii) the two-stage semi-randomized model proposed in this 
section above. 

In the two-stage model we take time t* = 2.0 for transition from a contamination randomized 
stage to the mean-field stage with random initial condition. The expected number of infectives in 
node 2 at time t = 2.0 is 100 which is large enough and at the same time much smaller than the 
node population. We take L = 10 4 for number of realization in the second stage to evaluate the 
moments. 

In the SIC approximation the randomized model presented in Table [2] comprises four con¬ 
sequently independent MCs. Processes 1 and 2 represent independent outbreak in node 1 (1271) . 
Processes 3 and 4 represent migration of its infectives to node 2 (l45l) . Processes 5 and 6 represent 
migration of host susceptives from node 2 to node 1 and their possible removal due to contam¬ 
ination; they are analogous to MC (1481) with N 2 substituted by S 2 to be valid for all the stages. 
Analogously processes 7 and 8 represent contamination of S 2 i in the first node and migration of 
appeared infectives J 2 i to their host node (l49l) . Finally processes 9-10 represent the outbreak in 
node 2 ; they are analogous to MC with an additional flux (l29l) with the accurate account of the 
number of succeptives S 2 at all the stages. 

The results of computation are presented in Figures |9] and [TO] Here the bold lines indicate the 
full randomized model, the thin solid lines indicate the SIC approximated model, the line with dots 
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indicate the two-stage semi-randomized model. Also the mean-field solution is presented as well 
and indicated by dashed lines. 




Figure 9: Comparison of different approximation for the mean value of infectives (left) and its 
standard deviation (right) for the populations Ni = N 2 = 10k. Bold line - full randomized model 
(Table 1), thin line - approximate randomized model (Table 2), line with dots - the proposed 
two-stage semi-randomized model, dashed line - the mean field solution. 




Figure 10: The same as in Figure [9]but for Ni = N 2 = 2k. 

Evidently, the proposed two-stage semi-randomized model gives a quite satisfactory approxi¬ 
mation for the first two moments of total number of infectives in node 2 if the population is 10k 
but only qualitative similarity for the smaller population. This justifies the used simplifications for 
rather moderate populated sites, but for the sites with population 2k and smaller a more sophisti¬ 
cated models should be developed. This can be a material of the consequent works. 


5 Discussion 

The randomized network SIR model coupled by randomized migration fluxes is described in terms 
of a Markov chain (MC). In the absence of infectives the pure migration model give reasonable 
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modelling of migration of individuals is described as a simple MC: being disturbed the system 
returns to a dynamic equilibrium exponentially fast that resembles a diffusion process in physics. 

In the mean-field (hydrodynamic) limit the MC converges to a non-standard network SIR 
model: the host and guest species are treated separately in the corresponding ODEs ©-(7]). 

Note that traditionally the coupling with other nodes is described by adding some transport 
terms (cf. Il6lfl2l0 


±S n = -/3 n S n I n + Xm n S, 

I n 3n^n^n ^n-^n 


Xmn^rn 


Simple analysis shows that pure migration in the such equation results in exponentially growing 
solutions 0. This instability is ignored as it can be hidden in the background of the outbreak and 
not be observable in certain epidemic model scenarios. 

The model proposed in 021 14] 


d q 
d t° n 

—i 

At ±n 


= —/3 n S n I n + \ S mn S '. 


V s w 

m A nm^ri 

Xmn^rn ~ Xnm^n 


s„ 


3 n S n I n oi n I n 


gives more stable pure migration but the simple analysis shows that in this model we obtain the 
fully mixed population in all the nodes 0 (e = 0.5 in our terms). Thus the dynamics of this 
model seems to be more realistic but nevertheless does not satisfy an intuitive interpretation of the 
equilibrium of the migration process. 

The both above models are characterized by only one parameter describing migration of a given 
species. This makes impossible to tune the model to obtain the realistic migration in the absence 
of the outbreak. 

In our earlier work [|4| a migration model is introduced without splitting species into host and 
guest with two migration parameters: £ and r describing migration process of a given species. But 
model proposed in 0 also does not give a completely satisfactory solution for pure migration in 
the case of different migration times: rf 2 and r 21 as shown in 0. 

Note that for some combination of the network parameters (epidemic and migration) effect 
of the more correct accounting of migration can be very small but it can become essential when 
parameters of the model vary. 

It is interesting to compare three different techniques for the model under consideration: a MC 
describing the number of individuals from all categories in both centers, its hydrodynamical limit 
in the form of a system of dynamic equations and a simple description of contamination stage at 
the node 2 as an isolated center with an inflow of infectives neglecting the backward migration. 
We developed an approach when the random evolution on the contamination stage either in its 
full or simplified form is coupled with the dynamical description on the stage of saturation. This 
makes the problem computationally feasible. Our intention is to apply this technique to a network 
of interacting population centers in the next paper. Note that the direct simulation of the network 
is far too expensive in terms of computational time compared with the integration of systems of 
dynamical equations. The computational time may be considerably reduced if the simulation is 
required during a relatively short contamination periods only. 

The spatial dynamics of multi-species epidemic models is widely discussed in the literature (see 
[ 15] ! I6} QT, 18 1 19] [2Qi 21 1 22] [23]] and the reference therein). For a purely deterministic model the 
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account of different dynamics of host and guest species on the epidemic speed was studied in 0. 
The simplifying assumptions make the analysis tractable by may not adequately reflect reality. It 
seems that a network of stochasticly interacting centres of the type discussed above may provide 
more realistic but still tractable setting. 

In the next paper we intend to derive the traveling wave characteristic equation (cf. 01,01,0) 
and explore analytically an numerically the dependence of the mean epidemic speed and its stan¬ 
dard deviation on the network parameters. 
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A Sketch of the proof to Proposition 1. 


First, the mean values of the Markov chain (MC) converges to the solution of initial value prob¬ 
lem ©-([9]) by the LLN. Phenomenological sketch is given Section [D and the rigorous proof is 
analogous to that presented in []24ll25 il. 

We also have to establish the convergence in probability. For definiteness consider J 2 (t) and 
apply Chebyshev’s inequality for any e > 0 


P 


hit) 

A 


m 


> e 


< var [I 2 (t)] 
- A 2 e 2 


Recall that A is the population scaling parameter (see Section©. So, it is enough to check that 


var [J 2 (t)] = 0(A), A —>■ oo. (58) 

It is demonstrated numerically in Section [3] (see Figure©. Actually we see that the normalized 
standard deviation decays as A r ~ l/2 , or equivalently, the non-normalized standard deviation grows 
as iV 1 / 2 (i.e., as A 1 / 2 ), that implies (l58l) . 

The rigorous argument runs as follows. Consider the processes in Table 1 which cause the 
change in number of infectives in node 2 and outline the fluxes of infectives. These processes are 


# 

Event 

Rate 

9,16 

J 2 —)■ / 2 + 1 

P 2 I 2 S 2 + /3 2 /i 2 S 2 + 

11,14 

J 2 —)■ J 2 — 1 

atih + ihh 


Here the flux terms are underlined, the remaining terms describe the MC based stochastic SIR 
model ©El. So the real flux can be defined as 


u (t) — /3 2 /i 2 S' 2 + <5( 2 / 2 1 — 721^2- 

We construct the process J 2 _ 


Event 

Rate 

I2 —■ 1 12 + 1 

fd 2 l 2 + v 

J 2 —)• I2 — 1 

07/2 


with the majorized constant flux 

d = P' 2 N\ + 8[ 2 N 2 > ift). 


( 60 ) 
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Remind that /3' 2 = /3 2 X 2 is a constant when A —* oo. 

For this process we have a randomized SI model (considered in |0) with the constant Poisson 
flux v. This problem is solved in Section 0 and it is shown that its variance grows as 0(A). 

Next, we establish the second order stochastic domination (see lf26l for details) of process / 2 (f) 
by J 2 (t). In fact the following inequality holds for all x, t > 0 (cf. 11261 ) 


rx rx 

/ P(/ 2 (t) > u)du < / P(/ 2 (f) > u)du. 

Jo Jo 

The second order stochastic domination means that for any convex function T(.) we have the 
inequality for all t > 0 

E[*(/ 2 (t))] <E[fl/(J 2 (f))], 

/ 2 [t) is the number of susceptible in the tractable model described by (l60l) . In our case 'f (A") = 
(X — EX) 2 . This implies the inequality var [J 2 (f)] < var J 2 (f) . □ 
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